/J/V52 ~/3 -^ 2 / 

' ' ' ^ 3 a phi 


Research Institute for Advanced Computer Science 

NASA Ames Research Center 


A Deterministic Particle Method for 
One-Dimensional Reaction-Diffusion 

Equations 


Michael Mascagni 


(NASA-CR-199755) A DETERMINISTIC 
PARTICLE METHOD FOR ONE-DIMENSIONAL 
REACTION-DIFFUSION EQUATIONS 
(Research Inst, for Advanced 
Computer Science) 17 p 


N96-15853 


Unci as 


G3/64 0083947 


RIACS Technical Report 95.23 


November, 1995 


A Deterministic Particle Method for 
One-Dimensional Reaction-Diffusion 

Equations 


Michael Mascagni 


The Research Institute of Advanced Computer Science is operated by Universities Space Research 
Association, The American City Building, Suite 212, Columbia, MD 21044, (410) 730-2656 


Work reported herein was supported by NASA via Contract NAS 2-13721 between NASA and the Universities 
Space Research Association (USRA). Work was performed at the Research Institute for Advanced Computer 
Science (RIACS), NASA Ames Research Center, Moffett Field, CA 94035-1000. 


A DETERMINISTIC PARTICLE METHOD FOR 
ONE-DIMENSIONAL REACTION-DIFFUSION EQUATIONS 


Michael Mascagni 

RIACS, NASA Ames Research Center 
and 

Center for Computing Sciences, I.D.A. 

Abstract. We derive a deterministic particle method for the solution of nonlinear reaction- 
diffusion equations in one spatial dimension. This deterministic method is an analog of a 
Monte Carlo method for the solution of these problems that has been previously investigated 
by the author. The deterministic method leads to the consideration of a system of ordinary 
differential equations for the positions of suitably defined particles. We then consider the 
time explicit and implicit methods for this system of ordinary differential equations and we 
study a Picard and Newton iteration for the solution of the implicit system. Next we solve 
numerically this system and study the discretization error both analytically and numerically. 
Numerical computation shows that this deterministic method is automatically adaptive to 
large gradients in the solution. 


1. Introduction. 

We are interested in numerical methods for the solution of nonlinear reaction-diffusion 
equations in one species and in one spatial dimension. These equations can be studied as 
initial- value problems and have the form: 

(1) u t = vu xx + /(u), t > 0, , u(x, 0) = uo(^), v > 0. 

Partial differential equations (PDEs) of this form are often found in multi-species systems 
that model phenomena from chemical combustion to excitability in neurons [10, 7, 11]. 

A qualitative feature of solutions to (1) is that they often have sharp gradients and 
traveling fronts. Because of this, standard finite-difference methods often require very fine 
grids to resolve the sharp gradients. With a traveling front solution the gradient is steepest 
in a small region that moves with the front. An optimal finite- difference method would 
incorporate moving refined grids into a numerical solution algorithm. 

An alternative is to consider a Monte Carlo method for the solution of (1). Such an 
approach yields a particle-based method that is automatically adapting to sharp gradients, 
[18, 19, 20]. The major drawback with a Monte Carlo approach is that convergence is 
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stochastic and can be shown to behave as 0(N~ 1 / 2 ), where N is the number of particles, 
[21]. An advantage of the Monte Carlo approach is that its particle-based formulation can 
serve as a starting point for the construction of a deterministic analog. This paper explores 
such a deterministic particle method for monotonic solutions to equation (1). 

The plan of the paper is as follows. In §2 we briefly review a one-dimensional Monte 
Carlo method, called the gradient random walk (GRW) method, for (1) and recall it’s 
major qualitative and quantitative properties. In §3 we derive the deterministic analog 
for monotonic solutions to equation (1). This leads to an interesting system of ordinary 
differential equations (ODEs) that has nonlinearities corresponding to the linear diffusion 
term in (1) and is linear where (1) is nonlinear! In §4 we discuss some common boundary 
conditions for reaction-diffusion equations and techniques for incorporating them into this 
numerical method. In §5 we consider the forward and backward Euler methods for the 
solution of the resulting ODEs. The backward Euler method leads to consideration of 
a nonlinear tridiagonal system. A Picard iteration is shown to solve this system; and a 
Newton’s method solution is also analyzed. In §6 we present some numerical results for 
the Nagumo equation, a scalar reaction-diffusion equation with monotonic solutions that 
models nerve conduction. These numerical results motivate an analysis of the discretization 
error. This error can be further broken down into an error in the wave speed, that can 
be analyzed analytically, and a static discretization error, that can be studied numerically. 
Finally, in §7 we conclude by summarizing our results and discussing the extension of this 
method to multi-species reaction-diffusion systems and to higher dimensions. 

2. Monte Carlo Algorithms. 

We first review the gradient random walk (GRW) method [1, 5, 19] for equation (1). We 
assume that is monotonic, i.e., u(x,t ) < u(y,t ) when x < y, and that u( — oo,t) = 

0,u(+oo,f) = 1. The strategy is to represent v = u x by diffusing particles. We choose to 
use a particle representation for the gradient since integration of the gradient to recover 
the solution is a variance reduction technique. Differentiating (1) with respect to x, 

(2) v t = vv xx + f'(u)v, i?(x,0) = u' Q (x), 

and discretizing v as a sum of 6-functions with strength mj, we have 

N 

(3) v(x,t) = Y^mjS{x - Xj(t)), 

j = i 

where Xj(t),j = 1, . . . ,N represent the location of the N particles. We use capital X to 
indicate that the positions are random variables. The density of the particles determines 
the value of v, and heuristically one thinks of mj as the “mass” of the jth particle. 

We recover u from v by u(x,t ) = f_ v(x',t): dx'\ 


N 

u(x,t ) = rrijH ( x — Xj(t )) . 

7=1 


( 4 ) 
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Thus u is represented as a step function with jumps of size nij at Xj. If all the particles 
have the same mass, 77?., the value of u at x is in times the number of particles that lie to 
the left of x. The boundary condition at —00 is automatically satisfied, and the condition 
at +00 is satisfied on average, with fluctuations [19]. This corresponds to conservation of 
mass. 

As described in [5], there is much less noise in the computed value of u than of v because 
all the particles contribute to the value at any point, x. Moreover, if the particles have 
equal mass, their density is large precisely where u has large gradients. 

Once the initial data is discretized, the GRW method evolves the particle positions and 
masses such that u satisfies equation (1). This is done by a fractional step iteration in 
which the diffusion term is modeled by a random walk and the reaction term is modeled 
as exponential growth- or decay of the particle masses. For each time step the sequence is: 

A. Gaussian Random Walk Step: 

Xj(1 -f- At ) — Xj(t) A o j 

where the Oj are independent iV(0, 2uAt) random variables 

B. Evaluate uj = u(Xj(t + A t)),j = 1, . . . , N using (4) 

C. Kill or Replicate Particles: with probability |/'(uj)| At 

i. Kill particle if /' < 0 

ii. Create a new particle at Xj if /' > 0 

Note that the only way the particles interact with each other, and hence the only way that 
the nonlinearity of the equations is manifest in the algorithm, is that the local value of u 
depends on the positions of all the particles. Not surprisingly, this is the step that is most 
difficult computationally. Naively, computing equation (4) for all the Xj' s requires 0(N 2 ) 
work, but if the particles are sorted, it requires only 0(N ) work plus 0(N log N ) for the 
sorting. 

The variant of this method due to Chorin and Ghoniem [1, 12, 5] is essentially the 
same, except that instead of killing and replicating particles, the masses are increased or 
decreased according to the ordinary differential equation 



In the mean and to 0(At ) the two procedures are equivalent. 

Both methods have been used to solve traveling-wave problems in one-dimension. The 
mechanism of wave propagation is transfer of mass from behind the front to ahead of it. 
With randomized particle death and replication, particles are killed off behind the wave 
front, where /'(•) is mostly negative, and replicated ahead of the front, where /'(•) is mostly 
positive. With deterministic particle growth and decay, there is a graded transmission of 
mass from neighbor to neighbor. 

3. The Deterministic Algorithm. 

Among both the different Monte Carlo algorithms we discussed in the previous section, 
the Chorin- Ghoniem variant already contains deterministic elements. However, it requires 
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using time dependent masses for the particles. What we seek is a totally deterministic 
algorithm where the masses of the particles remain constant and equal throughout the 
computation. This kind of algorithm has been applied to convection-diffusion equations 
before [3, 6, 13, 15, 16, 17], and we wish to extend these methods to reaction- diffusion 
problems. 

To begin our derivation, let us assume that we have a monotonic solution to (1), u(x,t ) 
with 0 < u(x,t ) < 1 for allt > 0 and u(—oo,t) = 0, u(-foo,<) = 1. Then we may define the 
functions x,(f) implicitly by saying they mark the locations of where the solution attains 
certain values as: 


(5) «(*.•(*)>*) = t> 0. 

Obviously we have that Xo(t) = — oo and x^(t) = +oo. We now seek equations of motion 
for the x,(-)’s given, that they satisfy (1). Let us take the total derivative with respect to 
time of (5): 

( 6 ) = 0 = u x (xi(t),t)xi(t ) + u t (xi (t),t). 

A rearrangement of (6) gives: 


(7) ii(t) = 


-1 


Ux(Xi(t), t) 


U t (Xi(t),t) = 


-1 


U x (Xi(t),t ) 


uu xx (xi(t),t ) + f(u(xi(t),i)) 


Equation (7) gives an expression for the velocity of the zth marker point, in terms 

of certain derivatives of u and values of /(«(•)). We have defined the values of u at the 
points £,•(<), so let us use those values to evaluate the terms above. Let us denote u(x,(t), t ) 
as U{. Using formally second order finite-difference approximations we first approximate: 


(8) 


-1 ^ Zj-n(t) - 

^x(®t(^)>^) ^t+1 Hi— l 




Next we approximate: 


u xx (xi(t),t) 


Uj + l -Uj Ui-Uj- 1 

X, + i(t)-X<(<) Xj(t)-Zj-i(t) 

X, + l(t) + X,(t) \ _ ( X,-(t) + X,-l(t) 


( 9 ) 


N (xi - (-1 

_*_( . 

2 ( ^H-i Xi— x 


•Ef-j-l % i i J 

" lr l 


[X. -Xj-! X, + i - X, J 


Finally it is obvious that 

( 10 ) 


((),*)) = / ( Jj ) , 
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so that the resulting equation becomes: 


( 11 ) 



This equation is very odd in the sense that its linear terms correspond to the nonlinear 
terms of the PDE (1) and visa versa. The first term on the right-hand side of (11) is highly 
nonlinear and corresponds to the linear diffusion term of (1), while the nonlinear f(u) term 
of the PDE results in the second (linear) term of (11). This is not, however, the complete 
system to consider. Usually equation (1) is part of an initial- value or initial boundary-value 
problem. Thus we must consider ways of incorporating the various boundary conditions 
associated with these problems into this curious formulation. 

4. Boundary Conditions. 

The system of ODEs in (11) is quite contrary with regard to linear and nonlinear terms. 

It is also a contrarian with regard to boundary conditions since the most natural problem 
to discretize is the pure initial-value problem with a monotonic solution. The pure initial- 
value problem on the entire line presents difficulties for finite-difference and finite-element 
methods. In these cases either explicit methods are used and only that subset of the line 
containing the dynamic parts of the solution is included in the computation, or an artificial 
boundary is used. With this deterministic particle method no such ad hoc approach is 
required. 

Another way that this discretization is contrary to finite-difference and finite-element r: ’ 

methods is that here we perform a uniform discretization of the dependent variable, u. In 
finite-difference and finite-element methods it is space, an independent variable, that is 
discretized to give a system of ODEs. While this traditional approach makes certain types 
of error analysis more convenient, discretization of the solution permits computational 
element, the particles, to concentrate automatically in areas where the solution is changing 
most rapidly. 

For the monotonic pure initial-value problem we have that xo(t ) = — oo and a;yv(f) = +oo 
at all times, and so x 0 (t ) = x;v(t) = 0. For i = 2, . . . , N — 2 we can just apply equation (11). 

The difficulty arises for i = 1 and i = N — 1 since equation (11) involves the points rco and 
x n in the terms involving /. These difficulties can be avoided if we go back to the derivation 
for these two equations and rederive them, being mindful to avoid the singularities. We 
may think of the ODEs in (11) as a result of combining the ODEs derived from the particle 
method discretization of the two equations (a) = uu xx and (b) Ut = f{u ) under the 

assumption of a monotonic solution. Discretization of (a) yields ii = u 

which is valid at i = 1 and i = N — 1, as the terms involving the infinite quantities are zero. 

To rederive (b) let us use the relevant part of equation (7). This means consideration of 
*.•(*) = ■ At * = 1 we may approximate u x (xi(t),t) » gffj- = JV(z2 1 _ gl) , while 

at n ~ N — 1 we have u x (x^-i(t),t) « n ^ Xn XN 2 y These along with the contributions 



c 


A DETERMINISTIC PARTICLE METHOD 


from (a) give the following system for the pure initial- value problem: 

x 0 = 0, 


X\ = —V 

(12) x, = v 

XN - 1 = v 

xn = 0 . 


X 2 - Xi 

1 


-N(x 2 -xi)f ( — j , 


l 


N (i 

(xi+i — X{— i ) / f ^ ) , * = 2, . 


• ,A r -2, 


2JV-1 — XAf-2 


- N (xN-i - XN-2) f 


'N - 1' 
N 


Similarly, Dirichlet boundary conditions are. easy to impose with this particle method. 
Consider first a time independent Dirichlet boundary condition: 

(13) u(0, t) = Uo, u(l,t) = U u UoKUl 

Here the solution is to rescale the problem so that Uo = 0 and U\ = 1 and use the conditions 
xo = 0 and xn = 1. The situation is a bit more complicated when the Dirichlet boundary 
conditions are time-dependent: 

(14) u(0,t) = U 0 (t), u(M) = tfi(<)> U 0 (t)<U 1 (t), t> 0. 


Since we have discretized u in our particle method, time- dependent boundary conditions 
require the creation and destruction of particles along the boundary. For example, if Uo(t) 
is increasing and in some time interval gains N~ 1 in value, the left-most particle’s position 
must be set to zero. This has the effect of increasing the value of the solution at zero by 
iV _1 . Conversely, if Uo(t) is decreasing and diminishes by TV -1 in value, a new particle is 
introduced at zero. At x = 1 we use a similar strategy, except that we remove a particle 
when Uo{t) decreases and insert a particle when Uo(t) increases. 

In many cases Neumann boundary conditions are appropriate. The time-independent 
version of these conditions 


(15) u x (0,t) = U o , u x (l,t) = Ui, C/ 0 > 0, Ui > 0. 

Since we are restricting our method to monotonic solutions, we must have ^^(O,^) > 0 
and u x (\,t) > 0. Let us first consider the leftmost boundary condition. Suppose that the 
left-most particle is at the location x\ > 0. An estimate of the derivative is thus 

ui — u 0 1 1 


(16) 


«z(0,t) 


= Uo — ► xi = 


Xl Nx 1 1 NU 0 

Similarly, at the right-hand side assume that x/v-i is the right-most particle in the unit 
interval, then the right-hand Neumann condition becomes xn-i = 1 — 

It should be clear that time-dependent Neumann boundary conditions are really no 
more trouble. Indeed, if we have the conditions: 


(17) u r (0,t) = Uo(t), u x (l,t) = U 1 (t), U 0 (t)> 0, Ui(t) > 0. 

They are enforced via: X\ (t) = Nt j 0 ^ and o:a'-i(0 = 1 — jrdjtj • ® ne important point 

is that we are constrained to keep these particles in the unit interval so that 1 > X\ (t) = 
and 0 < x'n-i (f) = 1 — ■> which means that N > y and N > j for all 

time. 
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5. Time Discretizations. 

Let us now consider some simple time discretizations of (11). We will consider only 
the forward- and backward- Euler methods here, which are first order methods. However, 
the analysis of the backward- Euler method will indicate the magnitude of the stiffness 
expected in these ODEs when solved by other implicit methods, [4], 

Let us consider the pure initial-values problem that results in the system of ODEs in 
(12). Let us first give this equation the vector notation: x(t) = F(x(t)), where x(f) = 
(a:o(t), Xi (t), . . . , .T/v-i(t), xs(t)) T and F(x(t)) is the right-hand side of (12). We now 
introduce a time discretization with time step At and Xj = x(jAt) so that the forward 
Euler method is merely: 


(IS) 


Xn+1 X„. 

A t 


= F(x n ) -> x n+1 = x n + A tF(x n ). 


With our particular F(-), this equation is known to converge to the exact solution as 
At — > 0 as O(At), [2]. In fact, it can be shown that if . r , — 7 t < J- this 

v n 1 1 min,[(r,+i — x,- )(!,— — 2 v 

difference equation is numerically stable. This fact is very much analogous to the stability 
of the forward- Euler method for a finite- difference solution to these same equations, [8, 14]. 

Because of the serious time-step restriction, we should consider the fully implicit backward- 
Euler method: 


(19) 


Xn-f-i 

At 


= F(x v + i ) -» x„ +1 - A tF(Xr, + i) = x„. 


This equation forms a nonlinear tridiagonal system. To see this let us consider the ?'th 
equation of equation (19): 


(20) z” +1 - A tv 


L x t x i-l x i 


N 




The second term corresponds to the linear diffusion term in (1), and it can be simplified 

by reduction to a common denominator. Let us define Jv, n+1 = - m — jrtrr, — jrnr so 

hi -i ; .i Jhf+i ^ ) 

that equation (20) becomes: 


(21) 


zr +1 *"-i 1 + < +1 < +1 + = *?, 


= -At\K? +1 + -f 
d n+l _ 1 + 2At/\" +1 , 


u 


" +1 = -At 


N 


J{T+1 — _/ ( — 
1 o J \ tv 


Solving (21) requires dealing with a very odd nonlinear system. Even if the nonlinear 
terms from diffusion are held constant (the 7f" +1 ’s), the linear system in (21) is not 
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guaranteed to be diagonally dominant. 1 However, a simple Picard iteration that will 
converge is possible if we rearrange (21) to allow us to work with a diagonally dominant 
system. Let us put the terms involving /(•) on the right-hand side as follows: 


( 22 ) 


LI 


n+1 ,_»»■+ 1 


X 


■-7 + D] 


n+l „n+l 


+ V\ 


n+l„n + l _ 
X !+l _ 



n + 1 


— X 


n+l\ 
-1 ) 1 


L? +1 = -A tK? +1 , 
D t n+1 = 1 + 2AL/v" +1 , 
L/" +1 = — At/v " +1 . 


This nonlinear system is diagonally dominant, and so the Picard iteration defined as follows 
will converge, [9]: 

A. Initialization: x old = x n 

B. Iterate until converged: 

i. Solve the following for x new : 

L° ld = -AtK° ld , 

Df d = 1 + 2AtK? ld , 

(23) Uf u = -AtK? d , 

j^old^new jjold ^new _|_ jjold ^.new _ 

+ At Y f {if} (*■+■ - ) ■ 


ii. Set x old = x new 

C. Form the solution to (22) as x n+1 = x new . 

Notice that the linear system in (23) is of the same form as one encounters in a finite- 

difference solution to (12). In fact, instead of the term L° ld = —A tK° ld we would have 

~Ax^ ~ 111 the finite- difference case. These terms are very much analogous. Note that the 

denominator of — A tKf ld = t~^u — «i is the product of two neighboring variable 

v x « * I i-iK I i+i _I i ) 

spatial grid spacings, just like in the finite-difference case. 

Now let us consider solving (21) with a Newton iteration. Recall that the Newton 
method iteration for the solution to the nonlinear system G(x) = 0 is defined implicitly 
by the equation J(x 0 id)(Xnew — x 0 id ) = — G(x 0 j<;). Here the Jacobian matrix, J, has 
elements defined by [ J ] = ^ i . For equation (21) we have 


(24) <3,(x) = x, — i/At 


|_ j — j X j -|- i X i 


+ Atj- ( x i+l j -x?, 


1 Diagonal dominance guarantees both nonsingularity and the numerical stability of Gaussian elimina- 
tion without pivoting, the Thomas algorithm. Without diagonal dominance we must pivot, which slows 
down the computation. 
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where the final, converged, solution gives x n +i. Thus J has the form: 

j = i + l 

J = i 

j = i- l 
otherwise. 

The matrix J is tridiagonal with a very curious form: it is the sum of a strictly di- 
agonally dominant symmetric matrix and a real antisymmetric matrix. The result is a 
tridiagonal matrix that, in general, cannot be thought of as diagonally dominant. How- 
ever, for most cases of interest this system will be diagonally dominant but not symmetric. 
Since we are adding and subtracting the same number from the off diagonals, the ma- 
trix will remain diagonally dominant if [J]n+i < 0. This requirement is equivalent to 
-j? > ^ 2 ~f { 77 )' Re P lace the term involving / with its maximum on the unit 

interval, f max = max* € [ 0|1 ] /(x), and divide by A tN 2 to give p v(T ;+ l~~ x,)p > JNfmax- 
Now the term [N( ' i ' an approximation to the square of the first derivative. Since 

the solution we seek is monotonic, this is always bounded away from zero by a constant, 
H~ 1 > 0, thus we rewrite this equation as N > . Since the right-hand side is con- 

stant, choosing N larger makes the inequality true; thus by choosing N sufficiently large, 
the system can be made diagonally dominant. 

6. Numerical Results. 

In this section we consider the numerical solution of Nagumo’s equation, [11], as part 
of a pure initial value problem: 

(26) u t = u xx + u(l — u)(u — a), t > 0, u(x, 0) = uo(x ), 1 > a > 0. 

In this pure initial-value problem setting, Naugmo’s equation has a stable traveling wave 
solution given by u(x,t) = ]+e _ (z 1 _ 9t)/ ^ i where the wave speed is 9 = \/2 [a — |). When 

a = \ we have 6 = 0, and the solution is a stable standing wave. These properties are 
easy to check numerically and make Nagumo’s equation a good candidate to explore the 
properties of a new numerical method for this class of problems. 

We begin by computing numerically the solution to Nagumo’s equation with our particle 
method using ito(x) = • This implies that u(x,t ) = ■ hr figure 1, 

below, we plot the computed solution at uniformly spaced times for the case where a = 0.75, 
with N = 64 and A t = 0.1. Notice that the particles move with the solution trackixrg the 
steep gradient. In addition, by the end of the computations every single particle has moved 
past the position of the rightmost particle in the initial configuration. This shows how this 
method is naturally adapting to the solution and is equitable in its distribution of particles. 

A qualitative observation that can be noticed by scrutinizing figure 1 is that the particles 
seem to be moving slightly away from the exact solution. This observation seems to indicate 
that actual wave speed of the particles may be slightly greater than 0, the wave speed of 
the exact solution. Fortunately, one calculation that we can carry out explicitly is the 


(25) 


i 


— v&t 1 WA( f ( j_\ 
"T 2 •' V N ) ’ 


(li + l-X ,) 2 
1 + 


i/A < 1 i/A < 

(x i + i x 1 ) ^ ' (x,-x ,_,) 2 


-i/A/ 


NAt 


(xi — x;_ 1) 2 

l 0, 


/(*). 
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Particle solution of Nagumo’s equation, dt = 0.1 , a = 0.75 


°t=100.0, exact 

t=1 00.0, approx x 
t=75.0, exact _ 
t=75.0, approx a 

t=50.0, exact 

t=50.0, approx a 

t=25.0, exact 

t=25.0, approx A 

t=0.0, exact 

t=0.0, approx D 



-5.000 0.(])00 5.000 10.000 15.000 20.000 25.000 30.000 35.000 40.000 45.000 50.000 

computation of the exact wave speed of the particle method solution to (26) to compare 

to the wave speed of the exact solution. Consider the wave speed of the midpoint of the 
travelling wave. Recall that the midpoint of the solution to (26) is the point, x m id(t), 
defined implicitly by u(x m id(t),t) = |. By inspection we see that x m id(t ) = Ot in the 
exact solution. As a check, we note that the exact velocity of the midpoint is also given 
by applying equation (7) to the midpoint: 


(27) 


•Emid(t) — 


-1 


u x(^mid(t)t t) 


llzx 


i-I 

2 


1 

2 ~~ a 


It is straightforward to calculate that u x (x m id(t),t ) = and u xx (x m id(t),t ) = 0. This 


gives x mid (t) = —4\/2 




= \/2 (a — ^) = 6, as we expected. 

Let us now calculate what the wave speed is of the midpoint when the exact solution 
is substituted into the discrete equations. We must therefore consider applying equation 
(11), with i — mid, v = 1, and /(u) = it(l — u){u — a), to u(x,t ) = 1+e _ (x \ gc)/v r 2 • 
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First we show that the diffusion term in (11) zero for Nagumo’s equation. Assume we 
have N points in our discretization, where N is even. Hence mid = y. We also have that 


mid±l 


~ 1 1 (,)- 9t )/y 2 , which can be solved to give x m id+i(t) = 6t — \/21n( 


N 


mirf+l 


1). Similarly, x m id-i(t) = Ot — \/2 ln ( — 1): so that the differences in the diffusion 
term are given by: 


•Emtrf+l X m id 


-y/2 


In 


N 


(28) 


•Emid X m id — i — \/2 


-V2\n{ 2 ^ 

mid 

N 


mid + 1 
N i 


-A-*( 


N 


In 


-V^ln 


mid 

N_ 

lid 


- 1 

- 1 ) - In 

- 1 


\mid 


N 

mid. — 1 


- 1 


- 1 


N 


mid — 1 


Now we compute: 

(29) (Xmid-\-\ S-mid) (%mid 3-mtd — 1 ) = \^2 In 


N 


mid+1 


- 1 


N 

mid 


N 


mid — 1 


N 

mid 


The numerator of the term inside the logarithm can be rewritten as: 
(30) 

N \/ N \ N 2 N 


mid + 1 


1 


mid — 1 


-1 = 


N 


(mid + l)(mid — 1) mid — 1 mid + 1 


+ 1 = 


N 2 — N [(mid + 1) + (mid — 1)] 


+ 1 = 


N 2 - N [2 (mid)\ 


(mid + 1 )(mid — 1) (mid + l)(rmd — 1) 

JV^-JV [2(f)] +1 = 1 


+ 1 = 


(mid + 1 )(mid — 1) 

And thus (x m id + 1 — x m id) = (x m id — ^mid-i), thus making the diffusion term in (11) 
vanish. For convenience, in the sequel we will refer to the common value of this difference 
as h. 

Thus we are left to consider the equation 


(31) 


• N, w /l\ N n ..f 1. 

•i'mirf(t) — ^(^mid+l x mid—l)f ( 2 ) l o ) ’ 


The term y2 h is the reciprocal of the centered difference approximation to the first deriv- 


ative. Indeed it is the case that 

( 32 ) 

Tv w( Xmid T u(x rn { d 

2 h = 2 h 

1 


h 2 

— Ui(lmifld) T g ^XXX ( Xmid,t ) + 0(h 4 ) 


— 2 / - 

2 Ux(x m id,t) + irU xxx (x m id,t) + 0(1l 4 ) U x (x m id,t ) 


h 'ILxxx 

6~ u x (x mid ,t) +U[n } 
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For the exact solution one can show that u xxx (x m jd,t ) = — • Recall that above we 

stated that u x (x m id(t,),t ) = so that = — f- Substituting this into (31) 

gives: 


(33) 


t 1 + 24 +0(A<) 


= e 


i+^+Offt 4 ) 


Equation (33) implies that the computed solution should have (1) a velocity that is 
slightly greater than the actual solution, and (2) this difference in velocity goes to zero 
with h 2 . Intuition suggests that ho c so one can think of the error in the velocity as 

going to zero with the square of the number of particles. We notice both these phenomena 
in our computations. 



a 

N 

0.125 

0.25 

0.375 

0.5 

0.625 

0.75 

0.875 

4 

-1.71e-01 

-1.33e-01 

-7.32e-02 

O.Oe-fOO 

7.32e-02 

1.33e-01 

1.71e-01 

8 

-1.04e-01 

-7.36e-02 

-3.77e-02 

0.0e+00 

3.77e-02 

7.36e-02 

1.04e-01 

16 

-4.70e-02 

-2.68e-02 

-1.18e-02 

0.0e+00 

1.18e-02 

2.68e-02 

4.70e-02 

32 

-1.91e-02 

-8.78e-03 

-3.23e-03 

0.0e+00 

3.23e-03 

8.78e-03 

1.91e-02 

64 

-7.62e-03 

-2.80e-03 

-8.23e-04 

0.0e+00 

8.23e-04 

2.80e-03 

7.62e-03 

128 

-3.06e-03 

-9.00e-04 

-1.96e-04 

0.0e+00 

1.96e-04 

9.00e-04 

3.06e-03 

256 

-1.25e-03 

-2.91e-04 

-4.35e-05 

0.0e+00 

4.35e-05 

2.91e-04 

1.25e-03 

512 

-5.14e-04 

-9.58e-05 

-8.72e-06 

0.0e+00 

8.72e-06 

9.58e-05 

5.14e-04 

1024 

-2.13e-04 

-3.19e-05 

-1.43e-06 

0.0e+00 

1.43e-06 

3.19e-05 

2.13e-04 

2048 

-8.88e-05 

-1.07e-05 

-1.15e-07 

0.0e+00 

1.15e-07 

1.07e-05 

8.88e-05 

4096 

-3.71e-05 

-3.68e-06 

4.73e-08 

0.0e+00 

-4.73e-08 

3.68e-06 

3.71e-05 


TABLE 1. The error in the wave speed ( 6 ) in the particle method 
solution to Nagumo’s equation. Errors for various values of a and N are 
presented. 

To confirm this calculation we have computed the velocity of x m fd(<) for various values 
of a, and hence various values of 6 , and with various (power-of-two) numbers of points. 
We have tabulated our results in table 1. A logarithmic regression actually gives the error 
in the wave speed as 0{N~ 2 ) instead of the expected 0(N~ 2 ). 

We have now computed the systematic error in the wave speed of our numerical method. 
The overall error in our numerical method is the sum of this error, which should be 0(t), 
and the error in the solution is due only to the spatial discretization. Nagumo’s equation 
allows us to study this second type of error directly by examining the numerical solution 
when a = \ and 0 = 0. In this case we expect that the numerical solution will deviate from 
the exact solution only, due to what is “discretization” error for this peculiar numerical 
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Figure 2. Mean square error of the particle method solution for vari- 
ous numbers of particles. The solution was computed with a = | giving 
a wave speed of 0 = 0. 

method as the solution is a standing wave. In fact, we see in table 1 that in this case 
there is no error in the computed wave speed for all values of h. In figure 2 we plot the 
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asymptotic mean square error in the particle method solution versus the number of points 
used. Notice that the slope of the curve in figure 2, a log-log plot, is ss —1.5, meaning that 
the the error behaves as 0(iV _ ? ). Thus the overall error of our method is 0(N~^). 

7. Conclusions. 

We have derived an adaptive particle method for the solution of one-dimensional reaction- 
diffusion equations with monotonic solutions. This method is the deterministic analog of 
a Monte Carlo method for these same problems. Like the Monte Carlo method, this de- 
terministic method uses a discretization of the solution to obtain a method that naturally 
tracks the solution by resolving the gradient. Thus this method is particularly suited to 
equations which have very steep gradients in the solution. We have also analyzed the 
difference equations associated with this discretization and have proposed a Picard and 
Newton iteration for the solution of the nonlinear system of equations arising from the 
fully implicit backward Euler time discretization. We have also decomposed the error in 
the method into the error in computing the wave speed and the spatial discretization er- 
ror. This former error can be asymptotically computed, while the later has been studied 
numerically. The overall error in the method appears to be 0(N~ ?). This is much better 
than the 0(N~ * ) error in the Monte Carlo method. 

Many questions remain. The most important question being why is the accuracy 
0(N~ 2 ) when we have chosen a spatial discretization based on second order accuracy? 
Perhaps the difference equations themselves can be refined to give a more accurate formu- 
lation. Or, the method of imposing the boundary conditions was only first order. This 
may have contaminated the overall accuracy of the solution. Is there a better formulation 
for the boundary conditions? It is also important to consider how to extend this method. 
Three ways that one would like to extend this method is to include (i) nonmonotonic but 
positive solutions, (ii) multi-species systems, and (iii) several spatial dimensions. 
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